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Abstract: We study the structure and energetics of Auat clusters by means of parameter- 
free density-functional calculations (A*' < 8), jellium calculations (A'' < 60), embedded- 

C^ . atom calculations (A'^ < 150), and parameterized density- functional calculations (N < 40) 

in combination with different methods for determining the structure of the lowest total 
energy. By comparing the results from the different approaches, effects due to geometric 
' O ' packing and those due to the electronic orbitals can be identified. Different descriptors 

C ' that highlight the results of the analysis are presented and used. 
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1 Introduction 
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("^ ■ Clusters of gold atoms have become the maybe mostly studied class of clusters [1], partly due to the 

VO i possibihty to apply them in electronic devices pi , nanomaterials 3 and catalysis ^ . Despite this 

^P ■ popularity only little consensus has been reached concerning the structure of these clusters. Many 

studies devoted to this issue use combinations of experimental and theoretical methods [51 [7]. 

From a theoretical point of view, gold clusters offer an additional challenge due to the impor- 
tance of relativistic effects, most notably the strong spin-orbit couplings. Thus, in parameter-free, 
^j^ • electronic-structure calculations, one has to use special, relativistic potentials. With those, the 

P3 I smallest gold clusters are found to be planar El El QB 1111 1121 Ej , whereas their exclusion leads 

O ■ to three-dimensional structures |14[ ITS] . 

However, for not-too-small gold clusters, it becomes increasingly difficult to apply parameter- 
free, electronic-structure methods in the calculation of the properties of the gold clusters, partly 
K^ \ because the computational demands scale with the size of the system to at least the third power, and 

; I ■ partly because the number of metastable structures grows very rapidly with cluster size. Thus, in 

^ I parameter-free studies one often has to make significant assumptions on the structure of the system, 
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as for instance is the case in the study of Haberlen et al. ^Hl • One of the greatest disadvantages of 
the first-principles methods is their incapabihty of optimizing large number of randomly generated 
initial structures and thereby determining the true global total-energy minimum. One example 
of this is provided by Auy and Aug for which the structure was predicted in 2000 by Hakkinen 
and Landman |1()| to be a planar structure with T)2h symmetry and a three-dimensional capped 
tetrahedron, respectively. Three years later, Hakkinen et al. |S] showed that the lowest-total- 
energy structure of Auy corresponds to a planar structure consisting of a rhombus, capped with 
an additional atom on three of its sides, and that the ground state of Au^ is the same rhombus, 
with its 4 sides capped. In 2005, the results for Aus were confirmed by Walker J21 and Remade et 
al. PSI , whereas both works found a planar capped hexagon to be the global minimum of neutral 
Au7. 

Approximate methods may provide a useful alternative to the parameter-free methods. They 
are computationally less demanding, thus allowing for a detailed search in structure space so that 
structures for clusters with well above 100 atoms can be predicted in an unbiased way. On the 
other hand, being approximate it is not obvious how reliable they are. It is the purpose of this 
work to address this issue. To this end, gold clusters provide an excellent playground, partly 
because of the large uncertainty concerning their structure in combination with the large amount 
of studies on these clusters, but also partly because for clusters both geometric packing effects and 
electronic shell effects may be responsible for the occurrence of certain particularly stable clusters 
(the so-called magic numbers) . The approximate methods often make different approximations on 
the relative importance of these two effects. 

One class of approximate methods is formed by the embedded-atom methods (EAM) that only 
indirectly includes electronic effects and, therefore, first of all (but not exclusively) put emphasis 
on packing effects. Both the EAM [171113111, the Sutton-Chen ^j, the Murrell-Mottram ,15] , 
and the many-body Gupta potential [2011^12312211^1131 (that all share the property of including 
electronic effects only very approximately) have all been applied in unbiased structure optimizations 
for gold clusters with up to 80 atoms. 

One of the, maybe, surprising outcomes of these studies is that the results depend very sensi- 
tively on the applied method, i.e., on the (more or less) approximate description of the interatomic 
interactions and on the method for structure optimization, see, e.g., [211 [23- One reason may be 
a subtle interplay between electronic and geometric effects, i.e., that the particularly stable struc- 
tures of Auat clusters are dictated partly by the closing of electronic shells and partly by geometric 
packing effects. Here, the various potentials give different relative importance to the two effects. 

The purpose of this contribution is to discuss general methods for calculating the properties of 
clusters, using gold clusters as the prototype. In parallel we shall also discuss the special properties 
of the gold clusters, specifically, with special emphasis on the issue above, i.e., whether electronic 
or packing effects are important in dictating the particularly stable clusters, and, moreover, how 
the different more or less accurate methods perform in calculating the properties of the clusters. 

2 Methods 

2.1 Total-energy methods 

The smallest gold clusters with up to eight atoms were treated with parameter-free electronic- 
structure calculations using the GaussianOS program package j28j. We performed density- functional 
(DFT) calculations using the generalized-gradient approximation (GGA) of Perdew, Burke, and 
Ernzerhof [2nilSn|- These calculations treat, in principle, all types of interactions, i.e., electronic 
and geometric effects, at an exact level. 

In addition we also performed self-consistent, electronic-structure calculations on spherical clus- 
ters where only the 11 (5d and 6s) valence electrons per Au atom were treated explicitly, whereas 
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all core electrons and the nuclei were smeared out to a uniform jcllium background |31[I32| . This 
model focuses essentially only on electronic effects. 

As an alternative we also considered the embedded-atom method (EAM) in the parameteriza- 
tion of Voter and Chen [^Sl 1^ OH] ■ According to this method, the total energy for a system of A'' 
atoms is written as 
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i.e., as a sum of atomic components, each being the sum of two terms. The first term is the 
energy that it costs to bring the atom of interest into the electron density provided by all other 
atoms, and the second term is a pair-potential term. Both terms are assumed depending only 
on the distances between the neighbouring atoms, and do therefore not include any directional 
dependence. Accordingly, the EAM emphasizes geometrical effects, whereas electronic effects are 
included only very indirectly. 

Furthermore, we used the density-functional tight-binding method (DFTB) as developed by 
Seifert and coworkers |36[l37j . With this method, the binding energy is written as the difference 
in the orbital energies of the compound minus those of the isolated atoms, i.e., as 



771 t 



(with m being an atom index and i an orbital index), augmented with pair potentials, 

/ ^ ^7771 ,7772 \ 1-^7771 -'^7772 | j V,^/ 

77715^7712 

(with Rm being the position of the mth atom). In calculating the orbital energies we need the 
Hamilton matrix elements {Xmini l-ff 1X7712772) ^-nd the overlap matrix elements (X777ini 1X7772772)- Here, 
X777n is the nth atomic orbital of the mth atom. The Hamilton operator contains the kinetic-energy 
operator as well as the potential. The latter is approximated as a superposition of the potentials 
of the isolated atoms, 

V{r)^Y.^U\r-Rm\), (4) 

777 

and subsequently we assume that the matrix element (Xn7ir7i 1^771X7772772) vanishes unless at least 
one of the atoms mi and m2 equals to. Finally, the pair potentials Umi,m2 ^^'^ obtained by requiring 
that the total-energy curves from parameter-free density-functional calculations on the diatomics 
are accurately reproduced. 

2.2 Structure determinations 

We used several different methods in determining the structures of the clusters. In the GaussianOS 
calculations the clusters were so small that it was possible to determine the structures of the lowest 
total energy simply through searching in the structure space. In the jellium calculations there is 
per construction no structure and the system has a spherical symmetry. 

In the EAM calculations we optimized the structure using our own Aufbau/Ahhau method 
[381 I39L 1^ . The method is based on simulating experimental conditions, where clusters grow by 
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adding atom by atom to a core. By repeating this process very many times and in parallel also 
removing atoms from larger clusters, we can identify the structures of the lowest total energy. 

Finally, in the DFTB calculations we used two different approaches. In one approach the 
structures of the EAM calculations were used as input for a local relaxation, i.e., only the near- 
est local-total-energy minimum was identified. In another set of calculations, we optimized the 
structures using the so-called genetic algorithms 01 02 ESI- Here, from a set of structures we 
generate new ones through cutting and pasting the original ones. Out of the total set of old and 
new clusters those with the lowest total energies are kept, and this process is repeated until the 
lowest total energy is unchanged for a large number of generations. 

3 Results 

The smallest possible cluster is the Au2 molecule. For this we show in Tabled the calculated bond 
length and binding energy from the different methods in comparison with experimental values. 
Notice that the DFTB method has been parameterized to reproduce results from parameter-free 
density-functional calculations on precisely the dimer and is, therefore, for the dimer accurate. 

Table 1: A comparison between the experimental and the calculated bond length and binding 
energy of the Au dimer obtained with the parameter-free density-functional calculations and the 

EAM method. 

Au2 Re, A Efc, eV 

DFT 2.55 2.22 

EAM 2.40 2.29 

EXP 2.47 2.29 

Next we show in Fig.Qlthe structures of Auat clusters with 4 < iV < 8 as obtained with the DFT 
calculations, the DFTB method, and the EAM method. The results for the trimer correspond to 
an obtuse triangle in the DFT calculations, and an equilateral triangle for the EAM and the DFTB 
methods. The figure clearly illustrates the aspects we have discussed above, i.e., the optimized 
structures of Au clusters result from a subtle interplay between geometric and electronic effects. 
Thus, in the DFT calculations all clusters are planar, whereas in the EAM calculations they are all 
three-dimensional (3D). Moreover, it turns out that also relativistic effects are important. Including 
all relativistic effects (i.e., also spin-orbit couplings) all gold clusters form planar structures at least 
up to iV = 13 (see Ref. |Sj). If the spin-orbit coupling is neglected, one obtains 3D global minima 
already at A^ = 4. 

The jellium model excludes packing effects and treats exclusively electronic-shell effects. With 
Ts being the electron-gas parameter of the system of interest (i.e., the radius of a sphere containing 
one electron) it is well-known ^^ that particularly stable clusters (magic numbers) occur for 
regularly spaced spherical clusters whose radius differ by 

AR = 0.603r, (5) 

for not too small clusters. For even smaller clusters one finds magic numbers for clusters containing 
2, 8, 18, 20, 34, 58, 92, 132, 138, 186, 254, 338, . . . electrons (see, e.g., 133]). 

In Fig. 121 we show results from the jellium calculations on Au clusters, where it is assumed 
that each atom contributes with 11 electrons. Thus, the above-mentioned magic numbers are 
not reached for the gold clusters. Nevertheless, the total energy per atom shows a very regular 
behaviour that actually can be related to Eq. Q . By comparing the total energy with the orbital 
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Figure 1: The optimized structures for Au^r, 4 < A'^ < 8 from (upper row) the DFT calculations, 
(middle row) the DFTB calculations, and (lower row) the EAM calculations. 
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Figure 2: Results from the jcUium calculations on Au^v clusters with 1 < A < 60. The three 
panels show (a) total energy per atom, (b) the stability function, and (c) the orbital energies (thin 
curves) and the Fermi energy (thick curve) as functions of A^. 



D.Alamanova, Y.Dong, H.u.Rebman, M.Springborg, V.G.Grigoryan 



energies (also shown) we see that the local minima correspond to structures where the Fermi level 
makes a jump, i.e., where new electronic shells are being filled. 

In order to identify the particularly stable clusters we introduce the stability function, 



A2^(iV) = EtotiN + 1) + EtotiN - 1) - 2^tot(A^) 



(6) 



where Etot{K) is the total energy of the Auk system. A2E{N) has local maxima when Auat is 
particularly stable, i.e., when Etot{N) is particularly low compared with Etot{N — 1) and _Etot(-^ + 
1). This function possesses a number of maxima, as seen in the figure, i.e., for A^ = 3, 8, 11, 16, 
21, 22, 26, 28, 30, 36, 37, 38, 47, 48, 50, 58, and 59. 
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Figure 3: Properties of Auat clusters from the EAM calculations. The four panels show (a) the 
stability function, (b) the radial distribution of atoms, (c) the similarity function, and (d) the 
shape-analysis parameters, respectively. Lengths and energies are given in A and eV, respectively. 
In (d) the upper rows show whether the clusters have an overall spherical shape (lowest row), an 
overall cigar-like shape (middle row), or an overall lens-like shape (upper row). 



We shall now compare this purely electronic description of the clusters with those obtained using 
the other approaches for Auat. Fig.^lshows various properties from the EAM calculations on Auat 
clusters with N up to 150. Here, we used our Aufbau/Abbau method in optimizing the structure. 
The stability function, Fig.jSJa), has pronounced peaks at A^ = 13, 54, and 146, where the structure 
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Figure 4: As Fig.|3 but from the DFTB calculations with the structures of the EAM calculations. 



corresponds to an icosahedron for the first two and a decahedron for the last (actually, for A^ = 55 
we do not find a icosahedron, whereas the structure of A^ = 54 is the 55-atoniic icosahedron without 
the central atom). 

Further information on the structure is obtained from the radial distances of the atoms, defined 
as follows. First, we define the center of the Au^r cluster. 



1 ^ ^ 



(7) 



and, subsequently, we define for each atom its radial distance 



n = \Ri - Rn 



(8) 



In Fig.|2Ib) we show the radial distances for all atoms and all cluster sizes. Each small line shows 
that at least one atom for the given value of A^ has exactly that radial distance. The figure shows 
that somewhere around A^ = 10 a second shell of atoms is being built up, with a central atom for 
A^ = 13. Aroimd A^ = 54, there are only few values of the radial distance, i.e., the clusters have a 
high symmetry. Around N — 75 we see that a third atomic shell is being formed. 
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Figure 5: As Fig. 13 but from the DFTB calculations with the genetic-algorithm optimization of 
the structures. 



We have earlier found [301 that it was useful to monitor the structural development of the 
isomer with the lowest total energy through the so-called similarity functions. We shall study how 
clusters grow and, in particular, if the cluster with N atoms can be derived from the one with 
A^ — 1 atoms simply by adding one extra atom. In order to quantify this relation we consider first 
the structure with the lowest total energy for the {N — l)-atom cluster. For this we calculate and 
sort all interatomic distances, di, i = 1,2,- ■■ , — ^— — -. Subsequently we consider each of the N 
fragments of the A^-cluster that can be obtained by removing one of the atoms and keeping the 
rest at their positions. For each of those we also calculate and sort all interatomic distances d^, 
and calculate, subsequently. 



N{N-l)/2 



N{N - 1) 



D E i'^^-^^)' 



1/2 



(9) 



Among the A^ different values of q we choose the smallest one, gmin, and calculate the similarity 
function 

S = ^—^- (10) 

1 + qmin/Ul 
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{ui = 1 A) which approaches 1 if the Auat cluster is very similar to the Auat-i cluster plus an 
extra atom. This function is shown in Fig. |3{c). We see that the structural development is very 
irregular over the whole range of N that we have considered here, with, however, some smaller 
intervals where S is relatively large, for instance for TV slightly above 20. 

Finally, we shall consider the overall shape of the clusters. As we showed in our earlier report 
on Ni clusters 0Dj, it is convenient to study the 3x3 matrix containing the elements 

1 ^ 

^st = ^y {Rn,s ~ Ro,s){Rn,t — Ro,t) (H) 



U 



= 1 



with s and t being x, y, and z. The three eigenvalues of this matrix, laa, can be used in separat- 
ing the clusters into being overall spherical (all eigenvalues are identical), more cigar-like shaped 
(one eigenvalue is large, the other two are small), or more lens-shaped (two large and one small 
eigenvalue). The average of the three eigenvalues, (laa), is a measure of the overall extension of 
the cluster. For a homogeneous sphere with N atoms, the eigenvalues scale like N^^^. Hence, we 
show in Fig. Ofd) quantities related to laa but scaled by TV^^/'^. In this figure we also mark the 
overall shape of the clusters through the upper points with the lowest row meaning spherical, the 
middle row meaning cigar-shaped, and the upper row meaning lens-shaped clusters. Some clusters 
with an overall spherical shape can be recognized, that, simultaneously, are clusters of particularly 
high stability according to Fig.|3fa). 

By comparing with the results of the jellium calculations, we see that the EAM method predicts 
a completely different set of particularly stable clusters. Moreover, compared with the results of the 
parameter- free density- functional calculations, the EAM method tends to produce more compact 
clusters. 

We shall now turn to the results obtained with the DFTB method. In Fig. 01 we show results 
similar to those of Fig. 13 but obtained by letting the clusters of the EAM calculations relax locally 
using the DFTB method, whereas we in Fig. [3 show the results for the DFTB calculations where 
the structure was optimized with the genetic algorithms. 

The structural information from all three sets of calculations is very similar. The radial dis- 
tances show in all cases that a second atomic shell is being constructed starting from slightly above 
10 atoms. First above iV = 40 a third atomic shell is formed [see Fig. |3Ic)]. Moreover, the sim- 
ilarity function shows that clusters in a narrow window just above A^ = 10 as well as in a wider 
window above A = 20 resemble each other, independently of the theoretical approach. There are, 
however, some differences in the overall shape, as seen in the shape-analysis parameters. In the 
EAM calculations several roughly spherical clusters with A^ < 40 are found, whereas the symmetry 
is broken when including electronic effects with the DFTB method, so that in these calculations 
no cluster is found to have an overall spherical structure. Moreover, the structures of the DFTB 
calculations are overall slightly more compact than those of the EAM calculations. 

On the other hand, the stability function of Fig. QJa) shows a much more irregular behaviour 
than the one of Fig.|2fa), although the structures are very similar. This emphasizes that electronic 
effects indeed are important. By comparing with Fig.Efa) we see that through further relaxations 
the stability function becomes less irregular, although its variation (over a range of roughly 12 
eV) is significantly larger than the range of the stability function in Fig.Ofa) (roughly 2 eV). It is 
interesting to observe that in the jcUium calculations (where there is per construction no structural 
relaxation) the stability function spans a range of roughly 35 eV, which is comparable to the range 
we see in Fig. EJa). This means that electronic effects are very important for the stability of 
the clusters, but also that through structural relaxations the role of the electronic effects can be 
somewhat reduced. 

Finally, we see in Fig. Q] that the DFTB calculations lead to essentially the same (type of) 
structures as the EAM calculations, with A^ = 4 being the only significant exception. Here, the 
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DFTB calculations predict a planar structure (as is the case for the relativistic DFT calculations), 
whereas for larger values of N, non-planar structures result in all our calculations except for the 
relativistic DFT calculations. 

4 Conclusions 

In this work we have discussed the interplay of electronic and packing effects in clusters. We have 
used Auat clusters as a prototype in order to illustrate the effects of different types of approxima- 
tions on the description of the interatomic interactions. We have, moreover, demonstrated how 
carefully chosen descriptors can be constructed that clearly grasp the essential outcomes of the 
calculations. Moreover, we hope also to have demonstrated the complexities related to unbiased, 
accurate calculations of the properties of not-too-small clusters. 

Our results indicate the existence of a tendency for simple potentials that do not directly 
include effects due to electronic orbitals (i.e., due to directed chemical bonds) to prefer closed 
packed structures. Including the electronic orbitals the structures may become less symmetric, and 
when increasing the accuracy of the treatment of the electronic orbitals and, simultaneously, the 
structure, the geometric arrangement of the atoms becomes less and less closed packed. Moreover, 
the electronic effects lead to a much stronger variation in stability as a function of size. 

We have in this presentation focused on Auat clusters. These clusters are among the most 
studied ones and are simultaneously very difficult to treat theoretically (as discussed in the intro- 
duction), which was the motivation for the present work. Nevertheless, we believe that most of 
our conclusions are valid also for other types of metal clusters, thus emphasizing the importance 
to perform different types of calculations with different types of (approximate) descriptions of the 
interatomic interactions before making finite conclusions about the properties of the clusters of a 
specific element. 
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